function [stt_o,eta_o,smooth_o,statedr_mat,etadr_mat]=kalman_drall(y,xdata,Z,Cobs,T,R,Qmat,A,addin,ndr,flag_robust)
% ==================================================================
% KALMAN_DRALL
% 
% This function takes as inputs the matrices of a state space model with no
% error in the observation equation and the Kalman Gain, and F^(-1) matrices from 
% a forward run of the Kalman Filter in order to simulate a NDR draws
% from the state disturbances and the states using a disturbance smoother. 
% It also returns the filetered and smoothed states. 
% It works both with ESTDGSE and SSM type models. 
% Only for NON_DIFFUSE state space with constant variances 
%
% The models is 
% y(t) = Z*s(t) 
% s(t) = T*s(t) + A*xdata(t) + R*n(t)  V( n(t) ) = Q=Qmat'*Qmat   
%
% Where y is [NZ  1], xdata(t) [NE  1]; s(t)[NS 1] and n(t) [NX x 1] 
% 
% Related files: 
% KALMAN_DR1.m          Generate a single draw 
% KALMAN_DR1_fast.m     Generate a single draw using the * variables in
% Durbin and Koopman, Biometrika 2003 
% Inputs 
% ======
% y     nobs x NZ stacked vector of observations, nobs is the total number of observations  
% 
% All matrices but Q(t) are assumed to be time invariant to simplify the computation 
% The correct size of input Qmat then depends on whether Q is time varyong or  not. 
% not time varying    Qmat is NX x NX 
%      time varying   Qmat is NX x NX x T where
% 
% ADDIN.
% FLAG_PZERO
% FLAG_SZERO 
% FLAG_ROBUST   Performs an additional smoother to ensure can recover
% observables 
%
% Alejandro Justiniano and Giorgio Primiceri (c) April 2005 
% =============================================================
if nargin < 10;error('Need define NDR, number of draws');end 
if nargin < 11; flag_robust=0;end; 
if flag_robust~=0 && flag_robust~=1;error('Flag robust must be 0 or 1');end 

% 1. Initial check for dimension
% -------------------------------
[nobs,nz]=size(y);
[ny,nx] = size(R);
ns=size(T,1);
Q=Qmat'*Qmat;
RQ=R*(Qmat'); 

% 2. Removal of Auxiliary data and means
% ydem
% ---------------------------------------
if ~isempty(xdata);
    %flag_xdata=1;
    if size(xdata,1)~=nobs;error('X and Y must have same number of rows');end
    if size(xdata,2)~=size(A,2);error('Cols A must match Cols of X');end;
    if size(A,1)~=size(T,1);error('Rows A must match dimension state space');end;
    ydem=y-xdata*(A');
else
    %flag_xdata=0;
    ydem=y;
end

if any(Cobs~=0)==1;
    ydem=ydem-repmat((Z*Cobs)',[nobs 1]);
end;

% 3. Initialization
% ==================
% I. Non-diffuse
% PSHAT P0|0
% SHAT  s0|0 
if addin.flag_pzero==1
    pzero=addin.pzero;
else
    pzero=lyapunov_symm(T,RQ*(RQ'));
    pzero=0.5*(pzero+pzero'); 
    if ~isempty(find(isnan(pzero))~=0)
        error('PSHAT contains Nan entries')
    end
end
if addin.flag_szero==1;
    szero=addin.szero;
    if length(szero)~=ns;error('SHAT must be NSx1');end
else
    szero=zeros(ns,1);
end

% ====================================================
% 4. Filter and Smooth YDEM 
% Resulting matrices have _o for orignal 
% =====================================================
Ztr=Z'; %clear Z;
Rtr=R'; %clear R;
[stt_o,eta_o,smooth_o]=kf_forbackfast_sub(ydem',szero,pzero); 
stt_o = stt_o'; 


% ============================================= 
% 6.    Simulate the innovations and the data
% w+
% y+
% s+    alpha+
% =============================================
etadr_mat=zeros(nobs,nx,ndr); 
statedr_mat =zeros(nobs,ns,ndr); 
iterset=100; 
jj=1; 
for jj=1:ndr;
    if rem(jj,iterset)==0; 
        dispaj(['Completed iteration ',num2str(jj)]); 
    end 
    wplus=randn(nobs,nx)*Qmat;
    % This is y+ and alpha+ in DK.
    % SZERO_PLUS: first state drawn from a Normal (shat,pshat)
    [yplus,stplus]=simssmodel(wplus,T,R,Z,szero,pzero);
    [stt_plus,eta_plus,smooth_plus]=kf_forbackfast_sub(yplus',szero,pzero);
    etadr_mat(:,:,jj)=(eta_o-eta_plus+wplus);
    statedr_mat(:,:,jj)=(smooth_o-smooth_plus+stplus);
end
% =========================================================================
% =========================================================================
% =========================================================================

    % =====================================================================
    % Soubroutine 
    % KF_FORBACKFAST_SUB 
    % Same as KF_DK but in subroutine mode 
    % STT and PTT DO NOT have the inital S0 and P0 in the first entry 
    % =====================================================================
    function [stt,etamat,smooth_st]=kf_forbackfast_sub(yy,a_zero,p_zero)
    %------------------------------------
    % a.FORWARD FILTER
    % -----------------------------------
    vstar=zeros(nz,nobs);
    finmat=zeros(nz,nz,nobs);
    kpartg=zeros(ns,nz,nobs);
    logLnc=zeros(nobs,1);
    %yfor=zeros(nz,nobs);
    % Matrices with one additional entry (initialization)
    stt=zeros(ns,nobs+1);
    ptt=zeros(ns,ns,nobs+1);
    stt(:,1)  =a_zero;
    ptt(:,:,1)=p_zero;
    ii=1;
    for ii=1:nobs;
        [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vstar(:,ii),finmat(:,:,ii),kpartg(:,:,ii)]...
            =kf_dk(yy(:,ii),Z,stt(:,ii),ptt(:,:,ii),T,RQ);
    end
    stt=stt(:,2:end);
    ptt=ptt(:,:,2:end);
    
    %------------------------------------
    % b.BACKWARD FILTER
    % -----------------------------------
    etamat   =zeros(nx,nobs);
    smooth_st=zeros(ns,nobs);
    
    rstar=zeros(ns,1);
    rmat=zeros(ny,nobs);
    
    [rstar,etamat(:,end)]=smoothdis(rstar,Q,Rtr,Ztr,finmat(:,:,end),zeros(ns),vstar(:,end));
    smooth_st(:,nobs)=stt(:,end)+ptt(:,:,end)*rstar;
    rmat(:,nobs)=rstar;
    
    ii=nobs; 
    for ii=(nobs-1):-1:1;
        [rstar,etamat(:,ii)]=smoothdis(rstar,Q,Rtr,...
            Ztr,finmat(:,:,ii),((T-T*kpartg(:,:,ii)*Z)'),vstar(:,ii));
        smooth_st(:,ii)=stt(:,ii)+ptt(:,:,ii)*rstar;
        rmat(:,ii)=rstar;
    end
    
    % =============================================================
    % c. Check if flag_robust==1
    % Perform another disturbance smoother and check can
    % recover the observables
    % ==============================================================
    if flag_robust==1
        smooth_st2=zeros(ns,nobs);
        smooth_st2(:,1)=T*a_zero+p_zero*rstar;
        ycheck=zeros(size(yy));
        ycheck(:,1)=Z*smooth_st2(:,1);
        for ii=2:nobs;
            smooth_st2(:,ii)=T*smooth_st2(:,ii-1) + R*(etamat(:,ii));
            ycheck(:,ii)=Z*smooth_st2(:,ii);
        end
        ydiscrep=max((abs(ycheck'-yy')));
        if max(ydiscrep > 1e-4);
            error('SMOOTH innovations do not recover observables');
        end
        
        % =============================================
        % Check also difference between smooth states
        if max(max(abs( smooth_st2-smooth_st ))) > 1e-4
            error('SMOOTH states do not match')
        end
    end
    
    etamat=etamat';
    smooth_st=smooth_st';
    
    end

end 